<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of demhmc2</title>
  <meta name="keywords" content="demhmc2">
  <meta name="description" content="DEMHMC2 Demonstrate Bayesian regression with Hybrid Monte Carlo sampling.">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../menu.html">Home</a> &gt;  <a href="#">ReBEL-0.2.7</a> &gt; <a href="#">netlab</a> &gt; demhmc2.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../menu.html"><img alt="<" border="0" src="../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="menu.html">Index for .\ReBEL-0.2.7\netlab&nbsp;<img alt=">" border="0" src="../../right.png"></a></td></tr></table>-->

<h1>demhmc2
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>DEMHMC2 Demonstrate Bayesian regression with Hybrid Monte Carlo sampling.</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>This is a script file. </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre class="comment">DEMHMC2 Demonstrate Bayesian regression with Hybrid Monte Carlo sampling.

    Description
    The problem consists of one input variable X and one target variable
    T with data generated by sampling X at equal intervals and then
    generating target data by computing SIN(2*PI*X) and adding Gaussian
    noise. The model is a 2-layer network with linear outputs, and the
    hybrid Monte Carlo algorithm (without persistence) is used to sample
    from the posterior distribution of the weights.  The graph shows the
    underlying function, 100 samples from the function given by the
    posterior distribution of the weights, and the average prediction
    (weighted by the posterior probabilities).

    See also
    <a href="demhmc3.html" class="code" title="">DEMHMC3</a>, <a href="hmc.html" class="code" title="function [samples, energies, diagn] = hmc(f, x, options, gradf, varargin)">HMC</a>, <a href="mlp.html" class="code" title="function net = mlp(nin, nhidden, nout, outfunc, prior, beta)">MLP</a>, <a href="mlperr.html" class="code" title="function [e, edata, eprior, mse] = mlperr(net, x, t)">MLPERR</a>, <a href="mlpgrad.html" class="code" title="function [g, gdata, gprior] = mlpgrad(net, x, t)">MLPGRAD</a></pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="hmc.html" class="code" title="function [samples, energies, diagn] = hmc(f, x, options, gradf, varargin)">hmc</a>	HMC	Hybrid Monte Carlo sampling.</li><li><a href="mlp.html" class="code" title="function net = mlp(nin, nhidden, nout, outfunc, prior, beta)">mlp</a>	MLP	Create a 2-layer feedforward network.</li><li><a href="mlpfwd.html" class="code" title="function [y, z, a] = mlpfwd(net, x)">mlpfwd</a>	MLPFWD	Forward propagation through 2-layer network.</li><li><a href="mlpinit.html" class="code" title="function net = mlpinit(net, prior)">mlpinit</a>	MLPINIT Initialise the weights in a 2-layer feedforward network.</li><li><a href="mlppak.html" class="code" title="function w = mlppak(net)">mlppak</a>	MLPPAK	Combines weights and biases into one weights vector.</li><li><a href="mlpunpak.html" class="code" title="function net = mlpunpak(net, w)">mlpunpak</a>	MLPUNPAK Separates weights vector into weight and bias matrices.</li></ul>
This function is called by:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="demnlab.html" class="code" title="function demnlab(action);">demnlab</a>	DEMNLAB A front-end Graphical User Interface to the demos</li></ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre>0001 <span class="comment">%DEMHMC2 Demonstrate Bayesian regression with Hybrid Monte Carlo sampling.</span>
0002 <span class="comment">%</span>
0003 <span class="comment">%    Description</span>
0004 <span class="comment">%    The problem consists of one input variable X and one target variable</span>
0005 <span class="comment">%    T with data generated by sampling X at equal intervals and then</span>
0006 <span class="comment">%    generating target data by computing SIN(2*PI*X) and adding Gaussian</span>
0007 <span class="comment">%    noise. The model is a 2-layer network with linear outputs, and the</span>
0008 <span class="comment">%    hybrid Monte Carlo algorithm (without persistence) is used to sample</span>
0009 <span class="comment">%    from the posterior distribution of the weights.  The graph shows the</span>
0010 <span class="comment">%    underlying function, 100 samples from the function given by the</span>
0011 <span class="comment">%    posterior distribution of the weights, and the average prediction</span>
0012 <span class="comment">%    (weighted by the posterior probabilities).</span>
0013 <span class="comment">%</span>
0014 <span class="comment">%    See also</span>
0015 <span class="comment">%    DEMHMC3, HMC, MLP, MLPERR, MLPGRAD</span>
0016 <span class="comment">%</span>
0017 
0018 <span class="comment">%    Copyright (c) Ian T Nabney (1996-2001)</span>
0019 
0020 
0021 <span class="comment">% Generate the matrix of inputs x and targets t.</span>
0022 ndata = 20;                     <span class="comment">% Number of data points.</span>
0023 noise = 0.1;                    <span class="comment">% Standard deviation of noise distribution.</span>
0024 nin = 1;                        <span class="comment">% Number of inputs.</span>
0025 nout = 1;                       <span class="comment">% Number of outputs.</span>
0026 
0027 seed = 42;                    <span class="comment">% Seed for random weight initialization.</span>
0028 randn(<span class="string">'state'</span>, seed);
0029 rand(<span class="string">'state'</span>, seed);
0030 
0031 x = 0.25 + 0.1*randn(ndata, nin);
0032 t = sin(2*pi*x) + noise*randn(size(x));
0033 
0034 clc
0035 disp(<span class="string">'This demonstration illustrates the use of the hybrid Monte Carlo'</span>)
0036 disp(<span class="string">'algorithm to sample from the posterior weight distribution of a'</span>)
0037 disp(<span class="string">'multi-layer perceptron.'</span>)
0038 disp(<span class="string">' '</span>)
0039 disp(<span class="string">'A regression problem is used, with the one-dimensional data drawn'</span>)
0040 disp(<span class="string">'from a noisy sine function.  The x values are sampled from a normal'</span>)
0041 disp(<span class="string">'distribution with mean 0.25 and variance 0.01.'</span>)
0042 disp(<span class="string">' '</span>)
0043 disp(<span class="string">'First we initialise the network.'</span>)
0044 disp(<span class="string">' '</span>)
0045 disp(<span class="string">'Press any key to continue.'</span>)
0046 pause
0047 
0048 <span class="comment">% Set up network parameters.</span>
0049 nhidden = 5;            <span class="comment">% Number of hidden units.</span>
0050 alpha = 0.001;                  <span class="comment">% Coefficient of weight-decay prior.</span>
0051 beta = 100.0;            <span class="comment">% Coefficient of data error.</span>
0052 
0053 <span class="comment">% Create and initialize network model.</span>
0054 <span class="comment">% Initialise weights reasonably close to 0</span>
0055 net = <a href="mlp.html" class="code" title="function net = mlp(nin, nhidden, nout, outfunc, prior, beta)">mlp</a>(nin, nhidden, nout, <span class="string">'linear'</span>, alpha, beta);
0056 net = <a href="mlpinit.html" class="code" title="function net = mlpinit(net, prior)">mlpinit</a>(net, 10);
0057 
0058 clc
0059 disp(<span class="string">'Next we take 100 samples from the posterior distribution.  The first'</span>)
0060 disp(<span class="string">'200 samples at the start of the chain are omitted.  As persistence'</span>)
0061 disp(<span class="string">'is not used, the momentum is randomised at each step.  100 iterations'</span>)
0062 disp(<span class="string">'are used at each step.  The new state is accepted if the threshold'</span>)
0063 disp(<span class="string">'value is greater than a random number between 0 and 1.'</span>)
0064 disp(<span class="string">' '</span>)
0065 disp(<span class="string">'Negative step numbers indicate samples discarded from the start of the'</span>)
0066 disp(<span class="string">'chain.'</span>)
0067 disp(<span class="string">' '</span>)
0068 disp(<span class="string">'Press any key to continue.'</span>)
0069 pause
0070 <span class="comment">% Set up vector of options for hybrid Monte Carlo.</span>
0071 nsamples = 100;            <span class="comment">% Number of retained samples.</span>
0072 
0073 options = foptions;             <span class="comment">% Default options vector.</span>
0074 options(1) = 1;            <span class="comment">% Switch on diagnostics.</span>
0075 options(7) = 100;        <span class="comment">% Number of steps in trajectory.</span>
0076 options(14) = nsamples;        <span class="comment">% Number of Monte Carlo samples returned.</span>
0077 options(15) = 200;        <span class="comment">% Number of samples omitted at start of chain.</span>
0078 options(18) = 0.002;        <span class="comment">% Step size.</span>
0079 
0080 w = <a href="mlppak.html" class="code" title="function w = mlppak(net)">mlppak</a>(net);
0081 <span class="comment">% Initialise HMC</span>
0082 <a href="hmc.html" class="code" title="function [samples, energies, diagn] = hmc(f, x, options, gradf, varargin)">hmc</a>(<span class="string">'state'</span>, 42);
0083 [samples, energies] = <a href="hmc.html" class="code" title="function [samples, energies, diagn] = hmc(f, x, options, gradf, varargin)">hmc</a>(<span class="string">'neterr'</span>, w, options, <span class="string">'netgrad'</span>, net, x, t);
0084 
0085 clc
0086 disp(<span class="string">'The plot shows the underlying noise free function, the 100 samples'</span>)
0087 disp(<span class="string">'produced from the MLP, and their average as a Monte Carlo estimate'</span>)
0088 disp(<span class="string">'of the true posterior average.'</span>)
0089 disp(<span class="string">' '</span>)
0090 disp(<span class="string">'Press any key to continue.'</span>)
0091 pause
0092 nplot = 300;
0093 plotvals = [0 : 1/(nplot - 1) : 1]';
0094 pred = zeros(size(plotvals));
0095 fh = figure;
0096 <span class="keyword">for</span> k = 1:nsamples
0097   w2 = samples(k,:);
0098   net2 = <a href="mlpunpak.html" class="code" title="function net = mlpunpak(net, w)">mlpunpak</a>(net, w2);
0099   y = <a href="mlpfwd.html" class="code" title="function [y, z, a] = mlpfwd(net, x)">mlpfwd</a>(net2, plotvals);
0100   <span class="comment">% Average sample predictions as Monte Carlo estimate of true integral</span>
0101   pred = pred + y;
0102   h4 = plot(plotvals, y, <span class="string">'-r'</span>, <span class="string">'LineWidth'</span>, 1);
0103   <span class="keyword">if</span> k == 1
0104     hold on
0105   <span class="keyword">end</span>
0106 <span class="keyword">end</span>
0107 pred = pred./nsamples;
0108 
0109 <span class="comment">% Plot data</span>
0110 h1 = plot(x, t, <span class="string">'ob'</span>, <span class="string">'LineWidth'</span>, 2, <span class="string">'MarkerFaceColor'</span>, <span class="string">'blue'</span>);
0111 axis([0 1 -3 3])
0112 
0113 <span class="comment">% Plot function</span>
0114 [fx, fy] = fplot(<span class="string">'sin(2*pi*x)'</span>, [0 1], <span class="string">'--g'</span>);
0115 h2 = plot(fx, fy, <span class="string">'--g'</span>, <span class="string">'LineWidth'</span>, 2);
0116 set(gca, <span class="string">'box'</span>, <span class="string">'on'</span>);
0117 
0118 <span class="comment">% Plot averaged prediction</span>
0119 h3 = plot(plotvals, pred, <span class="string">'-c'</span>, <span class="string">'LineWidth'</span>, 2);
0120 hold off
0121 
0122 lstrings = char(<span class="string">'Data'</span>, <span class="string">'Function'</span>, <span class="string">'Prediction'</span>, <span class="string">'Samples'</span>);
0123 legend([h1 h2 h3 h4], lstrings, 3);
0124 
0125 disp(<span class="string">'Note how the predictions become much further from the true function'</span>)
0126 disp(<span class="string">'away from the region of high data density.'</span>)
0127 disp(<span class="string">' '</span>)
0128 disp(<span class="string">'Press any key to exit.'</span>)
0129 pause
0130 close(fh);
0131 clear all;
0132</pre></div>
<hr><address>Generated on Tue 26-Sep-2006 10:36:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>